Lattice Boltzmann scheme for relativistic fluids 
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A Lattice Boltzmann formulation for relativistic fluids is presented and numerically verified 
through quantitative comparison with recent hydrodynamic simulations of relativistic shock-wave 
propagation in viscous quark-gluon plasmas. This formulation opens up the possibility of export- 
ing the main advantages of Lattice Boltzmann methods to the relativistic context, which seems 
particularly useful for the simulation of relativistic fluids in complicated geometries. 
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In the last decade, the Lattice Boltzmann (LB) method 
has attracted considerable interest as an alternative com- 
putational fluid dynamics method, based on the solu- 
tion of a minimal Boltzmann kinetic equation, rather 
than on the discretization of the equations of continuum 
fluid mechanics [1, 2]. To date, the overwhelming ma- 
jority of LB applications are directed towards classical, 
i.e. non-quantum and non-relativistic, fluids. However, 
while quantum versions of the LB equation have existed 
for more than a decade [11, to the best of our knowl- 
edge, an LB equation capable of handling relativistic flu- 
ids has not yet been proposed. In this Letter, we fill this 
gap and present an LB formulation for relativistic fluids. 
Our procedure is based on two simple and yet apparently 
unpursued observations, i) the kinetic formalism is nat- 
urally covariant /hyperbolic/conservative, ii) being based 
by construction on a finite- velocity scheme, existing lat- 
tice Boltzmann methods naturally feature relativistic-like 
equations of state, in the sense that the sound speed, c^, 
is a sizeable fraction of the speed of light c, i.e. the 
maximum velocity of mass transport {cg/c = with 
0.1 < < 1). Based on the above, one is led to propose 
that, upon choosing the lattice speed ci=5x/5t ~ c, the 
current LB mathematical framework should allow for rel- 
ativistic extensions, which is indeed the case as shown in 
this Letter. 

Our relativistic LB scheme (RLB) relies upon a 
moment-matching procedure similar to the one originally 
used to derive the earliest LB models for classical hydro- 
dynamics. That is, the local kinetic equilibria are ex- 
pressed as parametric polynomials of the relativistic fluid 
velocity jS = u/c^ with the parameters fixed by the con- 
dition of matching the analytic expression of the relevant 
relativistic moments, namely the number density, energy 
density and energy-momentum. As anticipated, the pos- 
sibility of a successful matching stems directly from the 
fact that, even in standard (non-relativistic) LB fluids, 
the sound speed Cg is of the same order of the speed of 




FIG. 1: Relativistic shock wave impacting on a sphere at 
1/3 1 = 0.5. Here the streamlines represent the velocity field, 
and the colors the pressure. The simulation was implemented 
using a grid size of 200 x 100 x 100 cells. 



light, typically Cg = c/V^, which is exactly the equation 
of state of relativistic fluids. As a result, = Ma/V^, 
so that is of the same order as the Mach number 
Ma = \u\/cs. Owing to this simple, and yet crucial prop- 
erty, it is therefore possible to tackle weakly relativistic 
problems much the same way as traditional LB handles 
classical low-Mach fluids. This spawns the exciting op- 
portunity of carrying the assets of LB over to the context 
of weakly relativistic fluids, such as the important case 
of quark-gluon plasmas generated bv recent experiments 
on heavy-ions and hadron jets 0, H 0, 0, S, i, [l^ • The 
RLB scheme is verified through quantitative comparison 
with recent one dimensional hydrodynamic simulations 
of relativistic shock-wave propagation in viscous quark- 
gluon plasmas [jjj. We can also apply our scheme to 
three dimensional geometries as shown in Fig. [H 

Being based on a second-order moment-matching pro- 
cedure, rather than on a high-order systematic expan- 
sion in P of the local relativistic equilibrium (Jiittner) 
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distribution [l2|, the present approach is in principle hm- 
ited to weakly relativist ic problems, with \P\ ~ 0.1. 
However, by introducing artificial faster-than-light par- 
ticles (numerical "tachyons"), the present RLB scheme 
is shown to produce quantitatively correct results up 
to \/3\ ~ 0.6, corresponding to Lorentz's factors 7 = 
> ^ ~ 1.4. Althoue^h still far from state-of-the-art 

numerical methods for relativistic hydrodynamics [13|, 
the RLB might nevertheless offer a fairly inexpensive al- 
ternative to more sophisticated methods at moderate val- 
ues of |/?|. In addition, since LB is recognizedly an excel- 
lent solver for fiows in complex geometries, like porous 
media, it is plausible to expect that the present RLB 
scheme may play a useful role for the simulation of rela- 
tivistic fluids in non-idealized geometries. 

To begin our model description, we focus on the rel- 
ativistic fluid equations associated with the conserva- 
tion of number of particles and momentum-energy. The 
energy-momentum tensor reads as follows 14, 15|: T^^ 



Pt]^^ + {e-\- P)u^u^ ^ with e the energy density and P the 
hydrostatic pressure. The velocity 4- vector is deflned by 
= (7, 7/?)^ , where P = u/c is the velocity of the fluid 

The tensor 



in units of the light speed and 7= > ^ 

77^^ denotes the Minkowski metric. Additionally, we de- 
flne the particle 4-flow, A^^ = (771,717/3)^, with n the 
number of particles per volume. Applying the conser- 
vation rule to energy and momentum, d^T^^ = 0, and 
to the 4-flow, d^N^ = 0, we obtain the hydrodynamic 
equations. 

Note that as opposed to a non-relativistic fluid, we 
have two scalar equations. To complete the set of equa- 
tions, we need to deflne a state equation that relates, at 
least, two of the three quantities: n, P and e. 

The above hydrodynamic equations can be derived as a 
macroscopic limit of the following relativistic Boltzmann- 
BGK equation [3, [H 
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where = {E{p)^pc) is the particle 4- momentum with 
E{p) the relativistic energy as function of the momen- 
tum modulus p=\p\^ f^^ a local relativistic equilibrum 
and r the relaxation time. Lattice Boltzmann theory for 
classical fluids shows that it may prove more convenient 
to solve fluid problems by numerically integrating the 
underlying kinetic equation rather than the macroscopic 
fluid equation themselves. The main condition for this to 
happen is that a sufliciently economic representation of 
the velocity space degrees of freedom be available. Fol- 
lowing upon consolidated experience with non-relativistic 
fluids, such a representation is indeed provided by dis- 
crete lattices, whereby the particle velocity (momentum) 
is constrained to a handful of constant discrete velocities, 
with suflicient symmetry to secure the fundamental con- 
servations of fluid flows, namely mass- momentum-energy 




FIG. 2: Set of discrete velocities for the relativistic lattice 
Boltzmann model. The highest speed is V^ci 



conservation and rotational invar iance. The main ad- 
vantages of the kinetic representation for classical fluids 
have been discussed at length [2], and they amount basi- 
cally to the fact that the information is transported along 
straight-streamlines (the discrete velocities are constant 
in space and time) rather than along space-time depen- 
dent trajectories generated by the flow itself, as it is case 
for hydrodynamic equations. Moreover, diffusive trans- 
port is not described by second-order spatial derivatives, 
but rather emerges as a collective property from adiabatic 
relaxation to local equilibria. This is crucial in securing a 
balance between flrst-order derivatives in both space and 
time, which is essential for relativistic equations. 

In order to reproduce the relativistic hydrodynamic 
equations, we propose an LB model with the D3Q19 cell 
conflguration, as shown in Fig. [2l We deflne two distribu- 
tion functions fi and gi for each velocity vector q, where 
the index i labels the discrete momenta within each cell. 
The hydrodynamic variables are calculated by impos- 



ing the following macroscopic constraints, = Xlilo 

(e + P)^2 _ p ^ ^18^^^^ (e + P)j^U = El^oa^C^■ 

From these equations, we have to extract the physical 
quantities e and P, where we have only flve equa- 

tions for six unknowns. The problem is closed by choos- 
ing an equation of state for ultra-relativistic fluids, e=3P. 

The distribution functions evolve according to the 
BGK Boltzmann evolution equation |l^ (full details in 
a future extended publication). 



U{x^m.t^5t)- U{x,t) =-^(/,-/;^) , (2) 



and. 



g^{x^m.t^5t)-g^{x,t) = -'-^{g, - g":^) , (3) 

where f^^ and g^^ are the equilibrium distribution func- 
tions. 

The equilibrium distribution functions that recover the 
relativistic fluid equations in the continuum limit, read 
as follows: 
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for i>0 and, 
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for OO and, 
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for the rest particles. Here, ci is the Umiting velocity of 
the lattice which relates the cell size and the time step 
Q = ||, and we have rescaled the velocity units such that 
the speed of light c=l. The weights for this set of discrete 
speeds are defined hy wq = 1/3 for the rest particles, 
Wi = 1/18 for the velocities |ci|=Q, and Wi = 1/36 for 
\ci\=V2ci. 

The choice of the state equation, e=3P, simplifies the 
equilibrium functions as follows. 
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for z>0 and. 
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P/Pq 0.70 



BAMPS: ri/s = 0.001 
BAMPS: ri/s = 0.05 
BAMPS: ri/s = 0.1 
LBS: ri/s = 0.001 
LBS: ri/s = 0.05 
LBS: ri/s = 0.1 




FIG. 3: Comparison between the BAMPS simulations [11] and 
the lattice Boltzmann results, for /3~0.2. Pressure (top) and 
velocity (bottom) of the fluid as function of the spatial coor- 
dinate z. 



for OO and. 
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for i=0. Then, the equations for the macroscopic vari- 



ables take the form: 



E1=o and 



4 2 ^ v^-L^ 



=0 9i^i' In our model, the shear 
viscosity is computed according to standard LB proce- 
dures as: 7^=|7^e(r — St/2)cf. 

To test the model we solve the Riemann problem in 
viscous gluon matter [11]. We use the state equation for 
ultra-relativistic fiuids e=3P, as before, and the rela- 
tion between energy density and particle number density, 
e=3nT, with T the temperature |l4j. The initial config- 
uration consists of two regions divided by a membrane 
located at z=0. Both regions have thermodynamically 
equilibrated matter with different constant pressure, Pq 
for z<0 and Pi for z>0. At t=0 the membrane is re- 
moved. 

We implement a one-dimensional simulation with an 
array of size 1x1x800. In this case, the 4- velocity is 
given by i^^=(7, 0, 0, 7/?)^. The velocity of the lattice is 
chosen q = 1.0, therefore the cell size Sx and time step 
6t are fixed to unity. This corresponds in IS units to 
^x=0.008fm and ^t=0.008fm/c. The viscosity is calcu- 
lated through 7^=|7^e(r — 1/2), and the entropy den- 
sity by the approximation s=4n — n In A, with the 



Initial condition - 
t = 0.64 fm/c-. 
t = 3.20 fm/c- 
t = 6.40 fm/c- 
t = 0.64 fm/c - 
t = 3.20 fm/c - 
t = 6.40 fm/c - 




FIG. 4: Time evolution of the shock wave for BAMPS 
simulations [11] and Lattice Boltzmann results, 
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FIG. 5: Velocity and pressure profile using numerical 
"tachyons" at time t=3.2fm/c, with f3^0.6 and ?7/s=0.001. 



gluon fugacity and the equilibrium particle density n^^ 
are given by, n^^= ^^^ with dG=lQ for gluons. Now, 
we can calculate the ratio between the viscosity and 
entropy density, 77/s, that is used as a parameter to 
characterize the shock- wave. The pressures were chosen 
Po=5.43GeVfm~^ and Pi=2.22GeVfm~^, corresponding 
to 7.9433x10"^ and 3.2567x10"^ lattice units, respec- 
tively. The initial temperature is To=350MeV, corre- 
sponding to To =0.0287 lattice units. With these parame- 
ters, the conversion between physical and numerical units 
for the energy, is lMeV=8.2x 10"^. 

Fig. [3] shows the results for different values of r]/s 
and the comparison with the BAMFSflTt (Boltzmann 
Approach of Multiparton Scattering) microscopic trans- 
port model simulations [11] at time 3.2fm/c. On the 
other hand, in Fig. [H we can see the evolution of 
the system for 7^/^=0.1 comparing the two numerical 
models. In both cases, we find an excellent agreement 
with BAMPS. To simulate fluid moving at higher speed, 
/3~0.6, we use numerical "tachyons" with q = 10. In- 
deed, from Eqs. [H and [5l it is seen that the positivity 

2 

condition /^^*^>0 implies q • u<^. Now, the pressure 
Pi is taken as 0.9532GeVfm~^ and we define two tem- 
peratures, To=0.0328 and Ti =0.0164, the first one for 
z < and the second one for 2: > 0. Fig. [5] shows 
the Shockwave for 7^/^=0.001 and the comparison with 
the BAMPS simulation [H], is again excellent. Our LB 
scheme easily extends to three dimensions, as is illus- 
trated in Figdl where we simulate the collision of a rel- 
ativistic shock wave with a fixed spherical obstacle. A 



typical 200 x 100 x 100 lattice-site simulations spanning 
1350 timesteps, takes about 1900 CPU seconds on a stan- 
dard PC. 

Summarizing, we have developed a Lattice Boltzmann 
formulation for (mildly) relativistic fluids. One of the 
major areas of application of non-relativistic LB schemes 
is flow through geometrically complex domains with in- 
ternal obstacles, such as porous media. It is therefore 
expected that the present RLB scheme may become use- 
ful for the simulation of relativistic fluids in complicated 
geometries. 
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